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Abstract 

In this paper we present an iterative method, inspired by the inverse iteration with shift technique of finite linear 
algebra, designed to find the eigenvalues and eigenfunctions of the Laplacian with homogeneous Dirichlet boundary 
condition for arbitrary bounded domains Q C R N . This method, which has a direct functional analysis approach, 
does not approximate the eigenvalues of the Laplacian as those of a finite linear operator. It is based on the uniform 
convergence away from nodal surfaces and can produce a simple and fast algorithm for computing the eigenvalues 
with minimal computational requirements, instead of using the ubiquitous Rayleigh quotient of finite linear algebra. 
Also, an alternative expression for the Rayleigh quotient in the associated infinite dimensional Sobolev space which 
avoids the integration of gradients is introduced and shown to be more efficient. The method can also be used in 
order to produce the spectral decomposition of any given function u 6 L 2 (Q). 

Keywords: Laplacian, eigenvalues, eigenfunctions, Fourier series, inverse iteration with shift, Rayleigh quotient. 

1 Introduction 

In Jl| we introduced an iterative method for computing the first eigenpair of the p-Laplacian operator A p u := 
div (jV«| p ~ 2 Vuj, p > 1, with homogeneous Dirichlet boundary condition in a bounded domain Q C R N , N ^ 1. 

The technique was inspired by the inverse power method or inverse iteration of finite linear algebra. 

In the present paper we concentrate in the special case p = 2, the Laplace operator A, which was superficially 
dealt with in HJ. Besides clarifying some of the arguments sketched in that paper for this case and providing some 
error estimates, our main purpose in this work is to show how inverse iteration with shift in the presence of uniform 
convergence can be used in order to obtain a fast and efficient method for computing the eigenvalues and eigen- 
functions of the Laplacian operator with homogeneous Dirichlet boundary condition for any bounded domain Cl. If 
the eigenvalues or at least good estimates for them are a priori known, the method can produce the corresponding 
eigenfunctions with great speed and accuracy. 

The technique can alternatively also be used as a fast process to obtain the spectral decomposition of any function 
u 6 L 2 (Q) (in other words, the Fourier series of u). 

We remark that the application of the method to the special case of the Laplacian operator is more natural since the 
Laplacian is a linear operator, L 2 (Q) is a Hilbert space and the inverse operator —A -1 is a self-adjoint and compact 
operator, therefore allowing the complete characterization of its spectral structure, as well as having the property that 
its eigenfunctions constitute a basis for L 2 (fi) (except for compactness, these properties are absent in the p-Laplacian 
when p 7^ 2). 

Our approach of the inverse iteration with shift is based on the following iterative process started by a given 
function u G L 2 (n) : 

<?o:=»and { = t CD 

T [ <pn+\ = ondfl 

'E-mail addresses: rodney@mat.ufmg.br (R. J. Biezuner), grey@mat.ufmg.br (G. Ercole), brenolg@ufmg.br (B. L. Giacchini) eder@iceb.ufop.br (E. 
Martins). 
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where o~ > is a previously fixed shift and A a := A + c\f is the corresponding shifted operator. 

The sequence {^«} is then handled in order to produce approximations for the pair (AJ,e£) where A^ denotes 
the eigenvalue of the Laplacian appearing in the spectral expansion of u which is closest to o~, and e£ denotes the 
eigenfunction obtained as the projection of u on the A^-eingenspace. 

Inverse iteration with shift is used in finite linear algebra in order to find the eigenvalues and eigenfunctions of a 
finite-dimensional linear operator. As an eigenvalue-finding procedure it is not as efficient as other methods, such as 
the QR algorithm. However, if very good estimates of the eigenvalues are known in advance, its rate of convergence 
to both eigenvalues and eigenfunctions can be very fast (see |2J, for instance). 

This approach can be naturally extended to self -adjoint compact linear operators in infinite-dimensional Hilbert 
spaces such as the Laplacian and those arising in Sturm-Liouville problems. In spite of this, we have not been able to 
find any reference in the literature to this approach being used in the Laplacian context. 

Since there is now a vast literature concerning the search for estimates for the eigenvalues of the Laplacian, as well 
as the gaps between eigenvalues (see, for instance, [3j |H[5l ED; although particularly useful in our context would be 
lower bounds for the difference between consecutive eigenvalues), these results can be used in connection with the 
inverse iteration with shift algorithm to find eigenfunctions of the Laplacian in arbitrary domains, as well as better 
approximations for its eigenvalues. It must be emphasized, however, that as with the finite linear method, the inverse 
iteration with shift method is not capable to find all the eigenfunctions associated to a non-simple eigenvalue. It can 
only find an eigenfunction of the associated eigenspace. In the generic sense most domains have Laplacian spectra 
consisting only of simple eigenvalues (see (ZHSl), although many domains of practical interest have eigenvalues with 
multiplicity greater than one (usually, domains which exhibit some type of symmetry, although not all of them). 

Algorithm 1 below is the simplest version of the inverse iteration with shift algorithm for computing one specific 
eigenvalue and corresponding eigenfunction of the Laplacian. 



Algorithm 1 Inverse Iteration with Shift for Laplacian Eigenvalue and Eigenfunction 

1: (p = U 

2: Set xq (point in Q outside nodal surfaces, randomly chosen) 

3: Set a (shift, usually eigenvalue estimate) 

4: Set m (number of iterations, depending on method used to solve Laplacian) 

5: for n = 0, 1,2, . . . , in do 

6: Solve — A^^n+i = tp n in O, (p„+\ = on 9Q 
7: end for 

8: return <f m +i/ ll'fm+illoo (L°° -normalized eigenfunction) 

9: return <p m (x )/(p m+1 {x ) + a (eigenvalue) 



In principle, the function u at the start of Algorithm 1 should be chosen so that it will have components in all 
eigenspaces of the Laplacian and a random choice would suffice. However, in practice, due to rounding errors 
any simple function can be used. In our numerical tests (see Section pk, we observed that even the unit constant 
function could be used in order to obtain all the eigenvalues, even in a domain where it does not have an infinite 
number of them in its spectral decomposition. 

In Line 6 of Algorithm 1 any PDE solver can be used. This allows one to choose the fastest solver available for a 
particular domain. In Line 9 the eigenvalue is computed according to the uniform convergence theory developed in 
Section]?] Since uniform convergence occurs away from nodal surfaces, a point xq £ O not in a nodal surface must be 
chosen; because nodal surfaces have zero N-dimensional measures, a random choice will suffice in the vast majority 
of cases, even taking into account that nodal surfaces change as the computed eigenfunction changes. 

In finite linear algebra, the approximation to the eigenvalue is usually computed using the Rayleigh quotient. In 
our context, the eigenvalue can also be computed via the Rayleigh quotient of the approximated eigenfunction: 

/ . v (Vfr,+i,Vfr,+i)2 llVjn+illj Jnl v ^+i| 2 dx 

><■ {<Pn + V = 7T T \ = — Tt" = J ■ ( 2 > 

(<p„+l,Cp„ + l) 2 \\(pn+l\\l Jn<Pn+l dx 

However, due to the high oscillatory nature of high frequency eigenfunctions, in order to accurately compute the 
integral of the (squared) gradient of eigenfunctions associated with these eigenvalues, a much finer grid needs to 
be used, which affects the efficiency of the method. Therefore, the Rayleigh quotient in its original form Q is not 
recommended for the computation of the eigenvalues of the Laplacian, unless one is prepared to incur the higher 
computational costs (see also further comments in Section ]7}. 

Nonetheless, an integration by parts argument produces the following alternative expression for the Rayleigh 
quotient that avoids computations involving gradients: 



Jci (pn<pn-ldx 

In \<l>n\ 2 dx 



2 



We obtain very satisfactory results by using this sequence in Line 9 of Algorithm 1 to compute eingenvalues in our 
numerical tests (see Section|6|. 

Another alternative way to compute the eigenvalues is given by the quotient 



Un+lh Jn<Pl+l dx 

when the shift a is chosen below the eigenvalue. This quotient also gives accurate approximations for the eigenvalues 
even using relatively coarse meshes; the computation of the integrals of the approximated eigenfunctions, instead of 
their gradients, does not appear to be significantly affected by the use of coarse grids (see Section [6|. 

We believe that, differently from what happens in finite dimensions, where much better and faster algorithms 
for finding the eigenvalues of linear operators (matrices), specially self-adjoint operators, are available, the inverse 
iteration with shift algorithm can be a very competitive method for finding the eigenvalues of the Laplacian. Typical 
algorithms for computing Laplacian eigenvalues involve the discretization of the Laplacian operator and the compu- 
tation of the eigenvalues of the resulting discretization matrix. However, since only a small number of the eigenvalues 
of the discretization matrix are good approximations to the Laplacian eigenvalues (the smaller ones), huge matrices 
are necessary in order to obtain a sufficiently good number of eigenvalues. And some problems, particularly those 
arising in the study of quantum billiards, demand the computation of a very large number of eigenvalues. Need- 
less to say, besides the requirements of memory, the size of the matrix makes it computationally costly to find its 
eigenvalues (see the classical [9| book, the review |3] and the more recent work |10J for details). 

The inverse iteration with shift method, that requires only typical relatively modest sizes for meshes in order 
to solve the Poisson equation with homogeneous Dirichlet boundary condition, can be very competitive in terms 
of memory allocation and processing time. This is true even when one considers that in order to find accurate 
approximations for the highest order eigenvalues and eigenfunctions sometimes one needs to refine the mesh, due 
to the increase of oscillations. 

Even if good estimates for the eigenvalues of a particular domain are not known in advance, a few iterations of 
inverse iteration with shift should be able to find good approximations to them, which can work as first estimates for 
the shift on a second run of the algorithm. The first choices for the shift might be concentrated around the numbers 
given by Weiyl's Law (see 1 11 J or Il2l ). 

The inverse iteration with shift method can also be used in order to find the spectral decomposition of any function 
u 6 L 2 (O), that is, in order to find its projections on the Laplacian eigenspaces. One has only to be careful to eliminate 
spurious projections, that is, projections which arise from rounding errors. This can be done through computing the 
Fourier coefficient associated to each eigenspace. If this coefficient becomes less than a specified very small tolerance, 
this projection can be safely discarded as arising from rounding errors. The two algorithms can be combined together 
in order to simultaneously find both the desired spectral decomposition of a given function defined on a domain and 
the spectrum of the Laplacian on it. The spectral decomposition algorithm is given below (Algorithm 2). Once again, 
Line 11 can be replaced by |3| or 

Algorithm 2 Spectral Decomposition 

1: (p = u 

2: Set xq (point in Q randomly chosen) 

3: for/c = 0,1,2,... do 

Set a k (shift) 
for n = 0, 1, 2, ... do 



Solve -Ao-^+j = <p* in Ci, <p* +1 = on 
Compute (u,<p^ +1 / \\<f$ l+1 \\ ) (Fourier coeficient) 



2/ 2 I 



12/2 

> tolerance then 



9: return <p* +1 / ||<^ +1 || (L 2 -normalized eigenfunction) 

10: return (u,<p k n+l / |^n+i| 2 ) (Fourier coefficient) 

11: return (p n (xo)/4>n+i{ x o) + &k (eigenvalue) 

break 

12: end if 

13: end for 

14: end for 



Although for the sake of simplicity all computations here are done for the Laplacian, the same algorithm can be 
used for similar elliptic operators. 
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This paper is organized as follows. The inverse iteration with shift sequence is defined in Section |2j where most 
of the notation used in this paper is also established. In Section 3 some well-known results concerning the Rayleigh 
quotient are recalled and proven for completeness. In Sections |3]and ijwe discuss L 2 and uniform convergence of the 
inverse iterated sequence, respectively. In the short Section|5]we make some considerations about the normalization 
process on each step of ([TJ which is useful for computational purposes, and in Section [6] we present the results of 
some numerical experiments made in simple domains. 

Finally, in Section[7]we discuss if the rate of convergence of the method could theoretically be improved through 
the use of the inverse iteration with shift given by the Rayleigh quotient, as it is standard practice in finite linear 
algebra. 

2 Definition of the inverse iteration with shift sequence 

Let O be a bounded domain in M N and H = {e k }\^ = ^ C Hq (Q) be an orthogonal (not necessarily normalized) 
basis for L 2 (fi) consisting of eigenf unctions of the Laplacian operator with homogeneous Dirichlet boundary condi- 
tion, that is, 

— Aej. = Aj-ej- in CI, 
ej. = on 3fi, 

where {^kjkLl * s me non-decreasing sequence of eigenvalues of the Laplacian, counting multiplicities: 

< Ai < A 2 < . . . (5) 

Let cr > and define the shifted operator 

Aa- :— A + crl. (6) 

It follows that e; c is also an eigenf unction of — A a corresponding to the eigenvalue A; c — a. 

Conversely, if A is an eigenvalue of — A a , then A = A^ — o~ for some k. Thus, the spectrum of the shifted operator 
equals the spectrum of the Laplacian operator shifted to the left by u, while the corresponding eigenspaces are the 
same. 

Given u £ L 2 (Q), let 

CO 

" = E K k e k 
k=l 

be the Fourier expansion of u, so that the Fourier coefficients are given by 

^ = (u,e k ) 2 = J n ue k dx 

Denote by A \ the least eigenvalue whose associated eigenspace is not orthogonal to u. That is, 

:= ^-h where k\ = min {k : cl\ 7^ 0} . 

In other words, A„ is the first eigenvalue such that u has a non-zero component in the corresponding eigenspace. 
Note that if r\ is the multiplicity of Ak then 

We will denote by e\ the orthogonal projection of u on the A \ -eigenspace, that is 

el ■= ocki% + H H ^h+n-i^+n-i = ^k- 

Thus, the Fourier expansion of u can be rewritten as 

" = E = 4 + E a ^k = 4 + E ^k- 

Proceeding in this way, denoting by e[, the orthogonal projection of u on the A[, -eigenspace which is the /^-eigenspace 
not orthogonal to u, the eigenfunction expansion of u can be written in terms of its non-zero components in the 
eigenspaces of the Laplacian as 

M . 

" = E e " ( 7 ) 

;'=i 
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where either M is a positive integer or, as in most cases, M = oo, and the corresponding sequence of eigenvalues 
| \[ t | is (strictly) increasing 

< x\ < a 2 < . . . 



As is well known from the theory of compact linear operators, if a does not belong to the spectrum of —A we 
have that (— A^)^ 1 : L 2 (Q) — > L 2 (fi) is a continuous, compact and invertible operator. Therefore, whenever a is 
not an eigenvalue of the Laplacian, we can define a sequence {(pn} ne ^ C Hq (fi) by inverse iteration setting 

cp = u and ( = *» in "' (8) 

For each w G L 2 (O) and a > not in the Laplacian spectrum consider the sequence {(p n } defined by inverse 
iteration in (jij. Since </> n +i = (— Au)^ 1 </>„, it follows from ^ that the eigenfunction expansion of (p„ is given by 

M j 

~ E Wi^' 
;'=! (a{, - a) 



Let A^J be the Laplacian eigenvalue appearing in the spectral decomposition of u which is closest to c, i.e., 

Mi 

Now let ef t denote the projection of u on the A^-eigenspace. Then 



(9) 



(u,<) 2 = «,<) 2 HK|| 2 (io) 



and we can write 



( 



M 

E 

A£-cr >|A£-crj 



M 

E 



(A{, " £T)« 



\ A; — (7 

V Ai— tr >|A5-£r| V ' 



or 



01! 



1 



n K + 'Pn), 



where 



M 

E 



Ai -O-j 



(11) 



(12) 



Throughout this paper, we will denote by X T U the Laplacian eigenvalue appearing in the spectral decomposition 
of u which is second closest to a, i.e., 



|A„ — crl = min 



X'u-cr 



(13) 



We will also denote by u n the component of u in the direction of ^ | , that is: 



u„ := < u, 



(pn \ (pn 



u ,, ,, dx 



\(pn\\ 2 /2 ll<M2 Il^"ll2 / 1 1 ^« 1 1 2 



(pit 



(14) 
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3 L 2 -convergence of inverse iteration with shift 

In this section we expose the L 2 -approach of the inverse iteration with shift. Firstly let us state some well-known 
results concerning the Rayleigh quotient TZ : H<J (fit) \ {0} — > K defined by 

For ease of consultation, we give short proofs of these results. 

Proposition 1 A function a 6 Hj (Ci) \ {0} is a critical point of TZ if and only if u is an eigenfunction of the Laplacian with 
homogeneous Dirichlet boundary condition and TZ (u) is the corresponding eigenvalue. 

Proof. Given v £ Hq (Ci), we have 

2 

TZ' {u)v = 2 K Vu ' Vj3 >2 ~ n (") ( u > v )i[ ■ 

||m|| 2 

Therefore, TZ' (u) = if and only if 

/ Vu-Vv = TZ(u) [ uv 
Jn Jci 

for all v e W Q U (Ci), that is, u is a weak solution of 

—Au = TZ (u) u infl, 
u — on dCl. 



Corollary 2 The Rayleigh quotient gives a quadratically accurate estimate for the Dirichlet Laplacian eigenvalues, that is, ifu is 
an eigenfunction of the Laplacian with homogeneous Dirichlet boundary condition with TZ (u) as the corresponding eigenvalue, 
then 

TZ(v)-TZ(u) =o(||i>-«|| 2 ) 

as v — y u in L 2 (Ci). 

Proof. If follows immediately from Taylor's formula, since TZ' (u) = 0. ■ 

Now we give the main result of this section. Before this, we note from 1 11 1 of Section [2] that 

<Pn ^ ± el + lpn 

\\<pn\\ 2 lk£ + <MI 2 ' 



(16) 



where the sign of the right-hand side will depend on whether the shift a is taken above or below the eigenvalue. 
With respect to n : if the shift is taken below the eigenvalue, the sign will always be positive, whereas if the shift is 
chosen above the eigenvalue the sign will be ( — !)"■ 



We note from 1 14 1 and ( 16 1 that 



u e a u +ip n \ e^ + ipn (17) 



K + WI2/2 K + <Mi2 

and also remark that 

Jnl«M dx 

This fact follows from an integration by parts after multiplying the equation —A a (p n = (p n -\ by <p n . From the compu- 



tational viewpoint, 1 18 1 has the advantage of avoiding the computation of the gradient V<p n as required in |[T5). 



Theorem 3 Let u G L 2 (Ci) 
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(i) We have 

Un\\ 2 < 

In particular, \p n — > in L 2 (O) with an exponential rate. 

(ii) There exists no £ N such that 



12 ■ 



"112 II 112 



for all n ^ no. In particular, 

with an exponential rate. 
(Hi) The following convergences hold: 

where u n is defined by pi] *, and 
with 



5r inL2{n) 

e u 112 



Vn ||2 



u„ — >• m L 2 (Q) 
fc(<pn) -> A£ 

ft(tf>„)-A£ = 



AJ-£T 



A? 



2« x 



Proof. We have from JT2b that 



IWI2 



M 

E 

Ui-<r >|Au— f| 



a4 - cr 



2h 



<: 



AST-cr 



A J 



2h 



E 

|A 7 „-(r|>|AJ-(r| 



<: 



A^-cr 



A J 



2/! 



12- 



Since 



< 1, 



it follows that 1 1 1 1 2 as n — » 00, which proves (i). 
Let n e N be such that 

1 

Wn\\ 2 < 2 " e "" for all n^n . 



Thus, if n ^ no it follows that 



^IKIl2=IKIl2-^IKIl2 



and 



^ IK + ^»|l2+ IWI2 ~ IWI2 = IK +^"112 

IKII2 + IK + 1M2 



kw + 'Pnlb KII2 



^2 



I e u + V'" 1 1 2 K 1 1 2 



e » (IKII2 K + V'" 1 1 2) "I" K II2 tyn 
(1/2)||<||1 

Kj|g IhMr+JKIjg H'M^ 

Hell 2 
IFm II2 



which proves (ii). 
Since 



\\<pn\ 



+I|l2 



IK + <MI 2 

\eZ + %+i\\z' 



pop follows from (z). 

The L 2 -convergence (21 1 follows from (27', (29> and (2CM. 
In order to prove ( 22 1 we firstly notice that 



lim 



(pn-1 <pn 



|| — 1 1| 2 ' ll^nlk/ 2 



-u 112 II K 112 / 2 



1 if A£ ^ CT 

-1 if A? < cr 



In fact, if A^ ^ c then 



lim 



<Pn-l <fa \ = Hm / g + gn-l < + ^ 
Il0n-l|| 2 ' 11^" II2 / 2 1 \IK + V" — 1 II2 ll e u tfn II2 1 2 



while if A„ < a then 



lim 



<Pn-l <Pn 



1 1 — 1 1 1 2 ' ll^ilk/ 2 



lim ((-if 1 /J + f*- 1 .. ,(-!)" e " + " / ' n 



Vn-ll 



2 / 2 



Thus, it follows from 1 18 1 that 



lim ft (<£„) = (cr+|A£-cr| 



1 if A^ t7 
-1 if < a 



A„. 



The convergence order | |23| follows from (21 1 and Corollary |2j since 



n (cp n ) - xi = n (u n ) - n «) = o (\\u n - ei\\ 



o 



2/i x 



4 Uniform convergence of inverse iteration with shift 

We begin this section by stating a L°°-estimate for an eigenfunction of the Laplacian in terms of its L 2 -norm. In 
the following, we denote by |0| the Lebesgue measure of Cl. 



Lemma 4 Let eeHj (O) be an eigenfunction of —A corresponding to the eigenvalue A. Then 

|H| M <4 N |n| 1/2 A N/2 || e || 2 . (24) 

Proof. It is shown in [13[, without any smoothness assumption on 3Q, that if e is an eigenfunction corresponding to 
a variational eigenvalue A of the homogeneous Dirichlet problem for the p-Laplacian then 

IkL^A^'iMi^. 

Choosing p = 2, pi) follows from Holder's inequality. ■ 

Estimates for eigenf unctions of the Laplacian with the exponent N /2 in the eigenvalue replaced by N/4 can be 
found in llT4llT5l[T6llT7l . See also [18, Remark 5.21] for more references. 

The following result refers to the nondecreasing sequence |5]l of eigenvalues of the Laplacian. 
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Lemma 5 If k > N/2, then 

OO -I j^^k 

where C is a positive constant which depends only on N and |C1| . 
Proof. It is well known (see OH, fM) that 

c 



A' > -/ 2/N 



where 

c «Mf! (26) 

N 47T 2 

and con is the volume of the N-dimensional unit ball. Hence, if / > N/2 we have 



i=i j=i 



2k -N 



In the next lemma we show that the convergence of a series formed by the eigenvalues A„ which appear in the 
spectral decomposition of a function u G L 2 (D) follows from the convergence of a series formed by all eigenvalues 
of the Laplacian. 



Lemma 6 Let k be chosen so that the series £ X is convergent. Then the series 

;'=l A / 



M (At) 



N/2 



A4 - «r 



N/2+fc+l 



zs a/so convergent. 

Proof. Assume that the expansion of u is not finite, i.e., M = oo (otherwise the result is trivial). Since 

OO OO -J 

;=i (a{,) ;=i j 



it suffices to show that 



(4) 



N/2 



N/2+jt+l 



* (Ai)' 



(27) 



for all sufficiently large j. 

As \ ! u — > oo, there exists jo such that 



A ; „ - o- 

N/2+k 



M, - o- 



N/2+k 



A ! u — o~ for all j ^ /o- Thus, if j is sufficiently large, we can write 

N/2+k 

<A i u -a 



(4 - tr) 



N/2+k 



M, - o- 



whence d27l follows. 



In order to prove the uniform convergence of the inverse iteration sequence {<pn} ne ?$, we return to 1 11 and write 

4 + 



<pn 
\\<Pn\\ t 



±- 



(28) 



As in (16) , the sign of the right-hand side will depend on whether the shift is taken above or below the eigenvalue 
and on n. 
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Lemma 7 The inequality 



(29) 



holds for all sufficiently large n, for some 6 > and a positive constant K = K (u,Ci, |A£ — cr\). In particular, xp n — > 
uniformly in CI with an exponential rate. 



Proof. From (12) and Lemma 1 we obtain 



M 

E 

\M, -(7 >|AS-£r| 



A ; H -<7 



N in 1 1/2 | 



|0| 



A4 

E 

Ui-cT- >|AS- 



A^ — cr 



a!, 



- (7 



(Ai) 



N/2 



But, taking = N/2 + £: + l, we have 



M 

E 



A£-<r 



Mi 



Ai 



N/2 



K-^ E 

A4-<r >|A£-£r] 



AS-<t 



a!, 



< |AS - *r|' 



A£-cr 



M 

E 





^N/2 


Ai- 


e 

- C 




N/2 


A{«- 




- t7 



AS-cr 



n-e m (x'u) 



N/2 



and by Lemma 3 the last series converges. Thus, 1 29 1 follows if we take 



K = 4 N |n| 1/2 || M || 2 K-^£ 



m (A{ 



N/2 



Ai 



A ; „-^ 



Theorem 8 Let u e L 2 (Q) . T/zen 
fz) X«ere exzsfs «o 6 N suc/z ffezf 

for all n ^ no- In particular, 

with an exponential rate, 
(ii) The following convergences hold 

and 



II <M 



uniformly in Cl 



ll^«+i l 



+ \K-'\ 

u„ — > e£ uniformly in Cl. 



(Hi) IfK C {x : e u u (x) ^ } is compact, then > A^ — a uniformly and with an exponential rate. 

<Pn+l 



(30) 
(31) 



Proof. Let no be such that H^Hoo ^ \ \\^u\\cx> f° r a ^ n ^ n 0- Then, as in the proof of (ii) of Theorem [3J we have for all 
n ^ no that 

1 „ 



10 



and 

K + M 

The remaining of (i) follows from Lemma |7| 
In order to prove ( 31 \ write 



\\<pn\ 



(pit 



\<Pn\\ 2 J ll^nL^n \\<pn 



I I'M 



<pn j 
u ,. ' — ax. 



Since 



and, from (i), 



lim 



(pu 



\\<Pn\ 



\\4>n\ 



lim 



K + *MI 

\K + lpn\ 



e u\\2 



ll^nL-/n \\<Pn 
uniformly in O, it follows from dTOk that 



<pn , e° u 
u ,, „ dx — > T , — -, 



(\K\ 


CO 


full / 


\ IK 


2 / 


' II^W II 00 & 



\p [r \\ Jo 



n \K\ 

.IT 



U T . — 7: — 



^dx 



-u 112 



Since from 111 I we have 



Nv4 

ll^n+l | 

and thus pO) also follows from Lemma [7] 
Now, let /C C C supp be compact so that 



= I A?, 



IK + ^»II C 



m :— min |e£| > 
K 



and fix uq S N such that 
Thus if n ^ ng we have on K, 
Therefore, the quotient 



m 

Un\L < 2 foralln > "0- 



K + $n\ ^ Kl - \lpn\ >nt - ll^nllao > 



H7 



0n 



makes sense on K for all sufficiently large n and again (Hi) follows from Lemma [7] since 



|AS-(r| 



1 



m 



e£ + ip n+1 



5 Normalization at each step 

In order to avoid numerical problems, such as overflow or underflow, it is usual to normalize the right-hand-side 
function in each inverse iteration. Although this procedure changes the sequence of iterates it maintains conver- 
gences. 

In fact, let v„ be defined by 



u 




Vq = and ||% || (32) 

on 3D. 
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where || ■ || may denote the L 2 -norm or the L°°-norm. Then, since a is not an eigenvalue of —A it is easy to verify that 

<Pn+\ - ||0«|| »n+l- (33) 
Hence, for example, if one uses || • || = || • || 2 then 



both uniformly and in L 2 . Note that this sequence is exactly the sequence {u n } defined by ( 14 and rewritten in terms 
of the sequence {v n } . 

Moreover, in view of (33) and (15) one also has 



n(v n ) = K{<p n ) ->■ K- (35) 
^ II ■ II = II ■ II oo then one has the following uniform convergence in each compact Kc{i: e u( x ) 7^ 0} : 

K ~ & f 1 if K ><r 



v n +i \K~<r\ I - 1 if A » <a - 

6 Numerical Tests 

In this section we present some numerical tests on the unit interval, unit disk and unit square. The inverse it- 
eration with shift was implemented starting from the unit constant function u = 1 on these domains. Eigenvalue 
approximations were computed by running our Algorithm 1 and taking (in the line 9) the following sequences con- 
sidered in this paper: 

<Pn , \ . \\<Pn\\2 
Pn :=-r^(xo)+£r, In ■■= || . II +a 
<Pn+l Wn+lh 

and 

TZ(4>n)=Cr^ " "2 ■ (36) 

ll<M 2 



This last sequence is an alternative form of the Rayleigh quotient evaluated at <p„ according to (18). As previously 
remarked, the numerical advantage of writing the Rayleigh quotient in this form is that one need not compute the 
gradient V(p„. 



Taking into account 135' we also compute eigenvalue approximations using the sequence lZ(v n ) where v n is 



defined by (32). Namely, we use the following alternative expression for this sequence: 

K(v n ) = a + -5— — - I LA . (37) 
IK|| 2 

Notwithstanding the (theoretical) equality TZ(v n ) = TZ((p n ), our tests indicate some remarkable numerical differences 
between them. 

The non-normalization of the function <p n on each step in Algorithm 1 tends to attribute to it smaller values at 
each iteration. Thus, as a numerical phenomenon, the quotient <p n l <p n +\ tends to assume the value 1, what makes the 
sequences }i n , j n and TZ((p„) converge to a + 1. 

This phenomenon may restrict the use of these sequences. However, handled with care by controlling the number 
of iterations and points in the grid, they can provide good approximations to the eigenvalues, as our numerical tests 
indicate. 

On the other hand, the sequence lZ(v n ) seems to be more robust, since we did not observe this phenomenon when 
using it. Another indication of its numerical stability is its tendency to better capture the correct eigenvalues (i. e. 
those belonging to the spectrum of the starting function), than the other sequences. 

The graphs of the eigenfunction approximations were constructed by using the sequence u n in (34) . In Algorithm 
2, these approximations can be viewed as the combination of lines 9 and 10 with replaced by (normalizing 
on each step). 

To show the efficiency of inverse iteration with shift we used neither the most efficient available method for 
solving the underlying differential equation nor a fine grid, but one of the most basic methods, finite differences, and 
a relatively coarse grid. Integrals were computed via the Simpson composite method. 
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Figure 1: First eight approximated eigenf unctions of the Laplacian on [0, 1] obtained from the inverse iteration with 
shift starting from the unit function. 

6.1 Eigenvalues and eigenfunctions for the unit interval [0, 1] 

In this case, (JSj becomes the following boundary value problem 

J ~<Pn+l ~ a< Pn+l = <Pn 

\ 4> n+1 (0) = 0„+i(i) =0. 

One can verify that Aj. = k 2 n 2 and that the function u = 1 does not have components corresponding to X\ for k even. 

We present in Table 1 exact and approximated eigenvalues of the Laplacian on the unit interval. The shift was set 
to the corresponding exact eigenvalue minus 0.1 (that is, u/ c := — 0.1), and a grid of only 101 nodes was used. As 
shown in the last column, the sequence lZ(v n ) seems to capture only the eigenvalues K a u that appear on the spectral 
decomposition of the function u = 1. In fact, we note that for even values of k the shift is closer to than A; c _j and 
even so the corresponding sequence lZ(v n ) converges to A^_^ which is the correct eigenvalue. 

In Figure 1 we present the graphs of the first eight approximated eigenfunctions. 



k 


A* 


Hw 


7io 




TZ(v 10 ) 


1 


9.8696 


9.8688 


9.8688 


9.8688 


9.8688 


2 


39.4784 


39.4654 


39.4654 


39.4654 


9.8691 


3 


88.8264 


88.7607 


88.7607 


88.7607 


88.7607 


4 


157.914 


157.706 


157.921 


157.706 


89.1623 


5 


246.740 


246.233 


247.047 


246.233 


246.233 


6 


355.306 


354.255 


356.157 


356.206 


252.178 


7 


483.611 


481.665 


485.356 


481.665 


481.665 


8 


631.655 


628.335 


634.773 


632.555 


514.785 



Table 1: Exact and approximated eigenvalues on the unit interval [0, 1] obtained from the inverse iteration with shift 
starting from the unit function. The shift <7) c := Aj- — 0.1 and a grid containing 101 nodes were used. 



Table 1 also exemplifies the numerical convergence to u + 1 of the non-normalized sequence 1Z((p n ), which hap- 
pened to the approximations of Ag and \g. For the first one, for example, the closest approximation achieved is 
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TZ{(pi) = 280.278. For n > 4 the quotient collapses to 1 and the result is spurious. In order to compute a correct 
approximation of this eigenvalue using this sequence, a finer grid should be used. 

In Table 2 we show the result of calculating the first 1,500 eigenvalues of the unit interval using the normalized 
sequence 1Z(v n ) with 30 iterations and a grid containing 10,001 nodes. The relative error between the computed 
eigenvalue lZ(v n ) and the exact eigenvalue A is defined by 



e(TZ(v n ),A) 



TZ{v n )-A 



A 



The shift used to make Table 2 is cry : = 0.99Aj-, which has an initial relative error of 1%. Such error is huge for great 
eigenvalues and the interval {cry, A/ c ) may contain many other eigenvalues. Thus, it is likely to happen that this shift 
makes the sequence converge to an eigenvalue A,- different from Aj-. With this in mind, we considered that TZ(v^q) 
correctly approximated an eigenvalue A r if |A,- — TZ(v^o)\ = min s |A S — TZ(v^q)\, and that TZ(v^o) converged to A r if 
their relative error is less than 10 -3 . 

Table 2 reveals that among the 1,500 shifts used, all of them approximated an eigenvalue with a relative error of 
order of magnitude 10 -3 , and that 1208 converged with a relative error less than 10~ 3 . 



e 


HT 3 


10" 4 


10" 5 


NT 6 


io- y 


^ 10" s 


N A 


292 


1004 


156 


39 


6 


3 



Table 2: Number of approximates with relative error of order e for ay := 0.99Aj-, grid containing 10,001 points and 
30 iterations. 

Among the 1208 converged approximations, 547 refer to eigenvalues Ay with k even. This shows that the sequence 
lZ(v n ) can converge to an eigenvalue that does not belong to the spectrum of the unit function. 

In order to better understand convergence of shifts with large initial relative error, we used a grid of 10,001 nodes 
and 30 iterations of the sequence TZ(v n ) with shifts cry := 0.5(Aj. + i + Aj.), where k runs from 1 to 100. The result is 
presented in Table 3. As we can see, despite the shift being located exactly between the eigenvalues, only one did not 
converge to an eigenvalue. 



e 10" 4 


10" b 


10" 6 


10~ 7 


^ 10" y 


N A 1 


66 


25 


4 


4 



Table 3: Number N\ of approximates with relative error of order e for cry := 0.5(Aj- + i + A/ c ), grid containing 10,001 
points and 30 iterations. 

Another numerical experiment using shifts with large initial error was done with randomly chosen shifts. We 
generated 100 random numbers (shifts) on the interval (0, A50) and used a grid of 10,001 nodes and 30 iterations. 
Table 4 shows the initial relative errors of the shifts, as well as the errors after 30 iterations of sequence lZ(v n ). As we 
can see, only one of these shifts did not converge to an eigenvalue, and most of them converged with a relative error 
of order 10 -5 . 



6 


^ 10-' 2 


10~ 3 


10" 4 


10" 5 


HT 6 


^ 10" 7 


No- 


38 


60 


2 











Na 





1 


2 


80 


16 


1 



Table 4: Numbers N t7 of shifts and N^ of approximates with relative errors of order e. Here a was randomly chosen 
on the interval (0, A50). A grid containing 10,001 points and 30 iterations were used. 



6.2 Radial eigenvalues and eigenfunctions for the unit disk 

We computed only the radial eigenfunctions for the unit disk Q = {xgR 2 : |x| <l}.In this case (p n = (pn(r) 
where r — \ x \ , and (ISl) becomes the Sturm-Liouville problem type 
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Km)' 



(T(f) n+1 =<p n , < r < 1 



Note that the function « = 1 has components in all radial eigenspaces. In fact, if e — e(r) denotes a radial eigenfunc- 
tion corresponding to an eigenvalue A > in CI then 

(re')' 

A r e = ± — '- = -\e < r < 1 



{ e'(0)=0 = e(l). 



Therefore, 



(w,l) 2 = / e(|x|)dx 

J |x| <1 

r)dS x dr = 2/r f 1 e(r)rdr = f\re'{r)Ydr = -^V(l) ^ 
JO A Jo A 



JLrl-r 



because of the uniqueness of the initial value problems for the ODE above at r = 1. 

We present in Table 5 the exact and approximated first eight radial eigenvalues for the Laplacian on the unit disk, 
calculated using the shift a k := \ k — 0.1 and a grid containing 201 nodes. 



k 






7io 




K(v w ) 


1 


5.7831 


5.7834 


5.7834 


5.7392 


5.7392 


2 


30.4713 


30.4698 


30.4698 


30.4396 


30.4396 


3 


74.887 


74.865 


74.865 


74.847 


74.847 


4 


139.040 


138.942 


138.942 


138.942 


138.942 


5 


222.932 


222.646 


223.018 


222.674 


222.674 


6 


326.563 


325.901 


327.025 


325.968 


326.563 


7 


449.934 


448.611 


451.057 


448.729 


448.729 


8 


593.043 


590.663 


595.223 


590.836 


590.836 



Table 5: The first eight Laplacian radial eigenvalues on the unit disk obtained from the inverse iteration with shift 
starting from the unit function. 

The graphs of the first eight approximated radial eigenfunctions of the Laplacian obtained by our Algorithm 2 
are displayed in Figure 2. 
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Figure 2: First eight Laplacian radial eigenfunctions on the unit disk obtained from the inverse iteration with shift. 
6.3 Unit square 

The eigenvalues of the unit square O = [0,1] x [0,1] are A n ,m = {n 2 + m 2 )n 2 and the corresponding L°°- 
normalized eigenfunctions are e, hm = sin(n7rx) sin(wz7rx). Hence it is easy to verify that the spectrum of the function 
u = 1 consists precisely of those eigenvalues \ n ,m for which both n and m are odd and that its first three eigen- 
values are X\ t \, A13 and A3 3. In Table 6 we present exact and approximated eigenvalues of the Laplacian on this 
domain. The shift was set cr^ : = Aj. — 0.1 and a grid containing 201 x 201 nodes was used. We can see again that 
the sequence TZ(v n ) tends to capture only the eigenvalues X a u that appear on the spectrum of the function w = 1, 
as shown in the last column. Note that the shift o\2 '■= M,2 ~ 01 is closer to A\ 2 but, however, the corresponding 
sequence lZ{v n ) approaches the correct eigenvalue Ai \ . The same behavior happens with the shifts 02,2 := A2,2 — 0.1 
and (72,3 := A23 — 0.1 since they are closer to A2,2 and A2,3, respectively, but the corresponding sequences TZ(v n ) ap- 
proach to A13. The graphs of the first three eigenfunctions in the spectrum of the unit function using Algorithm 2 are 
displayed in Figure 3. 



(n,m) 


Ajj,m 


m 


TlO 




TZ(v w ) 


(14) 


19.7392 


19.7388 


19.7388 


19.7388 


19.7388 


(1,2) 


49.3480 


49.3346 


49.3446 


49.3446 


19.8395 


(2,2) 


78.9568 


78.9504 


78.9504 


78.9504 


98.6732 


(1,3) 


98.6960 


98.6796 


98.6308 


98.6796 


98.6796 


(2,3) 


128.305 


128.285 


128.285 


128.285 


98.7042 


(3,3) 


177.653 


177.620 


177.620 


177.562 


177.620 



Table 6: Exact and approximated eigenvalues of the Laplacian on the unit square. 

In Table 7 we used the sequence fin to show the effect of refining the grid. The shift was chosen 0^3 ' — A3 3 — 0.1 
and 10 iterations were used. As expected, a finer grid provides a better approximation of the eigenvalue. 
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Grid 


Fio 


e(}t 10 , A 3/3 ) 


100 x 100 


177.5187 


7.6 x 10~ 4 


200 x 200 


177.6197 


1.9 x 10~ 4 


300 x 300 


177.6382 


8.3 x 10~ 5 


400 x 400 


177.6446 


4.7 x 10" 5 


500 x 500 


177.6476 


3.0 x 10~ 5 


1000 x 1000 


177.6516 


7.4 x 10~ 6 


2000 x 2000 


177.6526 


1.9 x 10~ 6 



Table 7: Approximated eigenvalues of the Laplacian on the unit square and relative errors for different grids. The 
exact eigenvalue is A3 3 = 177.6529. 




.A 



V v 



Figure 3: Graphs of the approximations for the first three projections of the function u = 1 on the unit square obtained 
of the inverse iteration with shift u„ m := X n , m — 0.1. e^ 1,1 (left), e^ 1,3 (center) and e^ 3,3 (right). 

7 Final comments 

In finite linear algebra, the iterative process itself is often used in order to generate increasingly better estimates 
for the eigenvalue at each iteration, meaning that the approximation obtained at any given iteration is used as the 
shift in the next iteration. It turns out that instead of using the estimates for the eigenvalue obtained in the process, 
the Rayleigh quotient of the estimates for the eigenvector obtained at each iteration give much better approximations 
for the eigenvalue. Indeed, if the eigenvalues of the operator or at least very good estimates of them are known in 
advance, inverse iteration with shift given by the Rayleigh quotient is the standard method for computing eigen- 
values due to its cubic rate of convergence (see |2|). It would be only natural to extend such ideas to the Laplacian, 
but we were not able to do it. Instead, our (admittedly preliminary) numerical tests, not shown in this paper, did 
not indicate convergence to the correct eigenvalues. As previously discussed, the Rayleigh quotient may not be a 
good way to approximate the eigenvalue of high frequency eigenfunctions unless the grid is much further refined, 
due to high oscillations, and the computational cost of using too fine grids can seriously limit the efficiency of the 
method. Further investigation is needed. So it remains an open problem to us if inverse iteration with shift given by 
the Rayleigh quotient is a method that can be successfully applied to the Laplacian. 

The method described in this paper uses a modification of the Rayleigh quotient that avoids the computation of 
gradients. Our numerical tests shown in Section 6 indicate a greater degree of convergence to the correct eigenvalues 
when this form is used. It seems to us that the non-necessity of calculating the gradients makes our algorithm more 
numerically stable. 

The only reference we could find where the Rayleigh quotient was used in computing the eigenvalues of the 
Laplacian, and only for polygonal domains, was the work [21J; however the Rayleigh quotient was only indirectly 
used there, as one component of another algorithm and in a very different way from the direct approach we follow 
here. 
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